library(tidyverse)
library(plotly)
library(sf)
library(tigris)
library(leaflet)
library(censusapi)
library(mapview)
library(png)
library(raster)

Sys.setenv(CENSUS_KEY="c8aa67e4086b4b5ce3a8717f59faa9a28f611dab")

Two separate sources of near-infrared (NIR) data

Load ecostress thermal data

The data comes from http://www.jnoble.org/ECOSTRESS/2020-08-03/baseMap1/index.html

ECOSTRESS Land Surface Temperature (LST) Silicon Valley, CA, USA – 2020.08.03, 9:42 pm PT Data: Glynn Hulley, ECOSTRESS team, NASA JPL

From layers.js in the .zip file of “Land parcels and buildings, 1:25,000 (OL web map, 230 MB), (.zip, 44 MB)” in http://www.jnoble.org/ECOSTRESS/index.html, I noticed the following details about the PNG:

projection: ‘EPSG:3857’,

imageExtent: [-13616154.835870, 4484784.913060, -13582758.988632, 4509845.252541]

So I looked up some techniques online, played around a bit, and was able to load in the PNG as an array, line it up with pixel coordinates in the right CRS, and then create a raster. Note I’m only using ecostress[ , , 2] which is only the “G” of the “RGB” data, which looks like it’s stored as fractions (we’ll need to figure out how to convert this back to thermal information, which should be something we can figure out from the literature on ecostress). There’s a flip() I need to do to get it in the right orientation – not exactly sure why but it looks like it’s in the right orientation at the end of this process. From here we can then do some work with the raster package, similar to how we intersect raster images with building footprints in this flood analysis: https://github.com/stanfordfuturebay/surge20/blob/master/SURF.Rmd#L80.

ecostress <- readPNG("G:/Shared drives/SFBI/CRNFO/ECOSTRESS_2020-08-03_SV_webmap1/ECOSTRESS_2020-08-03_SV_webmap1/layers/ECOSTRESSLSTC_2.png")

dat1 <- NULL

dat1$x <- seq(
  from = -13616154.835870,
  to = -13582758.988632,
  length.out = 932)

dat1$y <- seq(
  from = 4484784.913060,
  to = 4509845.252541,
  length.out = 700)

dat1$z <- ecostress[ , , 2] %>% t()

final <- raster(dat1) %>% 
  flip(direction = 'y')

crs(final) <- CRS('+init=EPSG:3857')

nfo_boundary <- 
  places("CA", cb = T, progress_bar = F) %>% 
  filter(NAME == "North Fair Oaks") %>% 
  st_transform(st_crs(final))

nfo_thermal_ecostress <- final %>% 
  crop(nfo_boundary)

mapview(nfo_thermal_ecostress)

Load SMC thermal data

nfo bounding box: 6063269, 1996269, 6072451, 2003697

maximum size image: 9182, 7428

https://gis.smcgov.org/image/rest/services/SanMateoCounty_Imagery2018/ImageServer/exportImage, manually query. request tif file. Can’t deliver full size, so my request was for half resolution, size 4591, 3714. Could ask for 4 quadrants separately and stitch together later if we want higher resolution.

library(raster)

smc_thermal <- raster("G:/Shared drives/SFBI/Data Library/SMC/NIR imagery 2018.tif")

nfo_boundary <- 
  places("CA", cb = T, progress_bar = F) %>% 
  filter(NAME == "North Fair Oaks") %>% 
  st_transform(st_crs(smc_thermal))

nfo_thermal_smc <- smc_thermal %>% 
  crop(nfo_boundary)

mapview(nfo_thermal_smc)

Even at half resolution, 2ft by 2ft. I’m not exactly sure what the units of this are, and if this is infact the “NIR” portion of the data (or just a camera image). What we have to go off of are https://gis.smcgov.org/image/rest/services/SanMateoCounty_Imagery2018/ImageServer, and literature out there. See https://www.earthdatascience.org/courses/earth-analytics/multispectral-remote-sensing-data/vegetation-indices-NDVI-in-R/, https://earthpy.readthedocs.io/en/latest/gallery_vignettes/plot_rgb.html for example.